println("Starting z_f with fixed theta.")

# make sure that we're using the right libraries
@everywhere using JuMP, KNITRO, FileIO, JLD2, Random

# do all prep work on every worker
@everywhere include("cenf.jl")

@everywhere CONST_DIM_ETA = 2

@everywhere begin

  root = "C:/Users/johannes.boehm/Dropbox/cenf/data/revision-2017"
  cenf.readData(root)

  # set seed, so that results can be replicated
  # my birthday
  rng = MersenneTwister(16091985)

  const CONST_MULTISTARTTRIES = 50

  ms_fvals=zeros(Float64, CONST_MULTISTARTTRIES,1);
  ms_θ=zeros(Float64, CONST_MULTISTARTTRIES,1);
  ms_β=zeros(Float64, CONST_MULTISTARTTRIES,2);
  ms_η=zeros(Float64, CONST_MULTISTARTTRIES,CONST_DIM_ETA);
  ms_γ=zeros(Float64, CONST_MULTISTARTTRIES,35*35);
  ms_logT=zeros(Float64, CONST_MULTISTARTTRIES,109*35);
  ms_logS=zeros(Float64, CONST_MULTISTARTTRIES,109*35);

  startβvec = Array{Array{Float64,1}}(undef,CONST_MULTISTARTTRIES);
  startηvec = Array{Array{Float64,1}}(undef,CONST_MULTISTARTTRIES);

  # generate starting values that satisfy constraint on beta
  for multiiter=1:CONST_MULTISTARTTRIES

    bDone = false
    while bDone==false

      startβvec[multiiter]=vec(rand(rng,Float64,1,2))
      # CONST_DIM_ETA
      startηvec[multiiter]=vec(1000.0 .*rand(rng,Float64,1,CONST_DIM_ETA));

      # this is the constraint on the β
      if (startβvec[multiiter][1]+startβvec[multiiter][2]<=1)
        bDone=true
      end
    end

  end

end


# do the search in parallel
# syntax:  pmap((a1,a2)->f(a1,a2),{"a","b","c"},{2,1,3})
@time outarray = pmap((startβ,startη)->cenf.solveWithFixedTheta(2.0,startβ,startη,1.0/35.0 .* ones(Float64,35,35),zeros(Float64,109,35),zeros(Float64,109,35)),startβvec,startηvec)


for i=1:CONST_MULTISTARTTRIES
    ms_fvals[i] = outarray[i]["fval"]
    ms_θ[i] = outarray[i]["θ"]
    ms_β[i,:] = outarray[i]["β"]
    ms_η[i,:] = outarray[i]["η"]
    ms_γ[i,:] = outarray[i]["γ"]
    ms_logT[i,:] = outarray[i]["logT"]
    ms_logS[i,:] = outarray[i]["logS"]
end
save("$(root)/output/ms_f_output_withfixtheta.jld2", "ms_fvals", ms_fvals, "ms_θ", ms_θ, "ms_β", ms_β, "ms_η", ms_η, "ms_γ", ms_γ, "ms_logT", ms_logT, "ms_logS", ms_logS)

@everywhere using JuMP, KNITRO, FileIO, JLD2, Random

@everywhere include("cenf.jl")

@everywhere begin

  cenf.readData(root)

  # second part:
  # search the neighbourhood of the best point for a better point
  f=load("$(root)/output/ms_f_output_withfixtheta.jld2")
  ms_fvals = f["ms_fvals"]
  ms_θ = f["ms_θ"]
  ms_β = f["ms_β"]
  ms_η = f["ms_η"]
  ms_γ = f["ms_γ"]
  ms_logT = f["ms_logT"]
  ms_logS = f["ms_logS"]

  idx_best = findfirst(isequal(maximum(ms_fvals)), vec(ms_fvals));
  bestθ = ms_θ[idx_best]
  bestβ = vec(ms_β[idx_best,:])
  bestη = vec(ms_η[idx_best,:])
  bestγ = reshape(ms_γ[idx_best,:],35,35)
  bestlogT = reshape(ms_logT[idx_best,:],109,35)
  bestlogS = reshape(ms_logS[idx_best,:],109,35)

  # set seed, so that results can be replicated
  rng = MersenneTwister(67890)

  const CONST_SECONDTRIES = 25

  ms_fvals=zeros(Float64, CONST_SECONDTRIES,1)
  #ms_θ=zeros(Float64, CONST_SECONDTRIES,1)
  ms_β=zeros(Float64, CONST_SECONDTRIES,2)
  ms_η=zeros(Float64, CONST_SECONDTRIES,CONST_DIM_ETA)
  ms_γ=zeros(Float64, CONST_SECONDTRIES,35*35)
  ms_logT=zeros(Float64, CONST_SECONDTRIES,109*35)
  ms_logS=zeros(Float64, CONST_SECONDTRIES,109*35)
  ms_r2 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varβ1 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varβ2 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varθ = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varη1 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varη2 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varη3 = zeros(Float64, CONST_SECONDTRIES,1)

  startβvec = Array{Array{Float64,1},1}(undef,CONST_SECONDTRIES)
  startηvec = Array{Array{Float64,1},1}(undef,CONST_SECONDTRIES)
  startγarr = Array{Array{Float64,2},1}(undef,CONST_SECONDTRIES)
  startlogTarr = Array{Array{Float64,2},1}(undef,CONST_SECONDTRIES)
  startlogSarr = Array{Array{Float64,2},1}(undef,CONST_SECONDTRIES)

  # generate starting values that satisfy constraint on beta
  for multiiter=1:CONST_SECONDTRIES

    bDone = false
    while bDone==false

      startβvec[multiiter]=bestβ + 0.2*(vec(rand(rng,Float64,1,2)) .- 0.5) .* bestβ
      startηvec[multiiter]=bestη + 0.2*(vec(rand(rng,Float64,1,CONST_DIM_ETA)) .- 0.5) .* bestη
      startγarr[multiiter]=bestγ + 0.2*(rand(rng,Float64,35,35) .- 0.5) .* bestγ
      startlogTarr[multiiter]=bestlogT + 0.2*(rand(rng,Float64,109,35) .- 0.5) .* bestlogT
      startlogSarr[multiiter]=bestlogS + 0.2*(rand(rng,Float64,109,35) .- 0.5) .* bestlogS

      # this is the constraint on the β
      if (startβvec[multiiter][1]+startβvec[multiiter][2]<=1)
        bDone=true
      end
    end


  end
end

# do the search in parallel
# syntax:  pmap((a1,a2)->f(a1,a2),{"a","b","c"},{2,1,3})
@time outarray = pmap((startβ,startη,startγ,startlogT,startlogS)->cenf.solveWithFixedTheta(2.0,startβ,startη,startγ,startlogT,startlogS),startβvec,startηvec,startγarr,startlogTarr,startlogSarr)

for i=1:CONST_SECONDTRIES
    ms_fvals[i] = outarray[i]["fval"]
    ms_θ[i] = outarray[i]["θ"]
    ms_β[i,:] = outarray[i]["β"]
    ms_η[i,:] = outarray[i]["η"]
    ms_γ[i,:] = outarray[i]["γ"]
    ms_logT[i,:] = outarray[i]["logT"]
    ms_logS[i,:] = outarray[i]["logS"]
    ms_r2[i] = outarray[i]["pseudor2"]
    try
      ms_varβ1[i] = outarray[i]["Var_β1"]
      ms_varβ2[i] = outarray[i]["Var_β2"]
    catch
      ms_varθ[i]= outarray[i]["Varθ"]
    end
    # CONST_DIM_ETA
    ms_varη1[i] = outarray[i]["Var_η1"]
    ms_varη2[i] = outarray[i]["Var_η2"]
    # ms_varη3[i] = outarray[i]["Var_η3"]
end
save("$(root)/output/ms_f_output_withfixtheta_refined.jld2", "ms_fvals", ms_fvals, "ms_θ", ms_θ, "ms_β", ms_β, "ms_η", ms_η, "ms_γ", ms_γ, "ms_logT", ms_logT, "ms_logS", ms_logS)


# VARIABLE THETA


println("Starting z_f, with variable theta.")

# Get more workers
# workers = addprocs(8)

# make sure that we're using the right libraries
@everywhere using JuMP, KNITRO, FileIO, JLD2, Random

# do all prep work on every worker
@everywhere include("cenf.jl")

@everywhere CONST_DIM_ETA = 2

@everywhere begin

  root = "C:/Users/johannes.boehm/Dropbox/cenf/data/revision-2017"
  cenf.readData(root)

  # set seed, so that results can be replicated
  # my birthday
  rng = MersenneTwister(16091985)

  const CONST_MULTISTARTTRIES = 50

  ms_fvals=zeros(Float64, CONST_MULTISTARTTRIES,1);
  ms_θ=zeros(Float64, CONST_MULTISTARTTRIES,1);
  ms_β=zeros(Float64, CONST_MULTISTARTTRIES,2);
  ms_η=zeros(Float64, CONST_MULTISTARTTRIES,CONST_DIM_ETA);
  ms_γ=zeros(Float64, CONST_MULTISTARTTRIES,35*35);
  ms_logT=zeros(Float64, CONST_MULTISTARTTRIES,109*35);
  ms_logS=zeros(Float64, CONST_MULTISTARTTRIES,109*35);

  startθvec = Array{Array{Float64,1}}(undef,CONST_MULTISTARTTRIES);
  startηvec = Array{Array{Float64,1}}(undef,CONST_MULTISTARTTRIES);

  # generate starting values that satisfy constraint on beta
  for multiiter=1:CONST_MULTISTARTTRIES

    startθvec[multiiter]=vec(10.0 * rand(rng,Float64,1,1));
    # CONST_DIM_ETA
    startηvec[multiiter]=vec(1000.0 .*rand(rng,Float64,1,CONST_DIM_ETA));

  end

end

# do the search in parallel
# syntax:  pmap((a1,a2)->f(a1,a2),{"a","b","c"},{2,1,3})
@time outarray = pmap((startθ,startη)->cenf.solveWithθ_f(startθ,startη),startθvec,startηvec)


for i=1:CONST_MULTISTARTTRIES
    ms_fvals[i] = outarray[i]["fval"]
    ms_θ[i] = outarray[i]["θ"]
    ms_β[i,:] = outarray[i]["β"]
    ms_η[i,:] = outarray[i]["η"]
    ms_γ[i,:] = outarray[i]["γ"]
    ms_logT[i,:] = outarray[i]["logT"]
    ms_logS[i,:] = outarray[i]["logS"]
end
save("$(root)/output/ms_f_output_withvartheta.jld2", "ms_fvals", ms_fvals, "ms_θ", ms_θ, "ms_β", ms_β, "ms_η", ms_η, "ms_γ", ms_γ, "ms_logT", ms_logT, "ms_logS", ms_logS)


@everywhere using JuMP, KNITRO, FileIO, JLD2, Random

@everywhere include("cenf.jl")

@everywhere begin

  cenf.readData(root)

  # second part:
  # search the neighbourhood of the best point for a better point
  f=load("$(root)/output/ms_f_output_withvartheta.jld2")
  ms_fvals = f["ms_fvals"]
  ms_θ = f["ms_θ"]
  ms_β = f["ms_β"]
  ms_η = f["ms_η"]
  ms_γ = f["ms_γ"]
  ms_logT = f["ms_logT"]
  ms_logS = f["ms_logS"]

  idx_best = findfirst(isequal(maximum(ms_fvals)), vec(ms_fvals));
  bestθ = ms_θ[idx_best]
  bestβ = vec(ms_β[idx_best,:])
  bestη = vec(ms_η[idx_best,:])
  bestγ = reshape(ms_γ[idx_best,:],35,35)
  bestlogT = reshape(ms_logT[idx_best,:],109,35)
  bestlogS = reshape(ms_logS[idx_best,:],109,35)

  # set seed, so that results can be replicated
  # my brother's birthday
  rng = MersenneTwister(12091983)

  const CONST_SECONDTRIES = 25

  ms_fvals=zeros(Float64, CONST_SECONDTRIES,1)
  ms_θ=zeros(Float64, CONST_SECONDTRIES,1)
  ms_β=zeros(Float64, CONST_SECONDTRIES,2)
  ms_η=zeros(Float64, CONST_SECONDTRIES,CONST_DIM_ETA)
  ms_γ=zeros(Float64, CONST_SECONDTRIES,35*35)
  ms_logT=zeros(Float64, CONST_SECONDTRIES,109*35)
  ms_logS=zeros(Float64, CONST_SECONDTRIES,109*35)
  ms_r2 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varβ1 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varβ2 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varθ = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varη1 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varη2 = zeros(Float64, CONST_SECONDTRIES,1)
  ms_varη3 = zeros(Float64, CONST_SECONDTRIES,1)

    startθvec = Array{Array{Float64,1},1}(undef,CONST_SECONDTRIES)
    startηvec = Array{Array{Float64,1},1}(undef,CONST_SECONDTRIES)
    startγarr = Array{Array{Float64,2},1}(undef,CONST_SECONDTRIES)
    startlogTarr = Array{Array{Float64,2},1}(undef,CONST_SECONDTRIES)
    startlogSarr = Array{Array{Float64,2},1}(undef,CONST_SECONDTRIES)

  # generate starting values that satisfy constraint on beta
  for multiiter=1:CONST_SECONDTRIES

      startθvec[multiiter]=bestθ .+ 0.2*(vec(rand(rng,Float64,1,1)) .- 0.5) .* bestθ
      startηvec[multiiter]=bestη .+ 0.2*(vec(rand(rng,Float64,1,CONST_DIM_ETA)) .- 0.5) .* bestη
      startγarr[multiiter]=bestγ .+ 0.2*(rand(rng,Float64,35,35) .- 0.5) .* bestγ
      startlogTarr[multiiter]=bestlogT .+ 0.2*(rand(rng,Float64,109,35) .- 0.5) .* bestlogT
      startlogSarr[multiiter]=bestlogS .+ 0.2*(rand(rng,Float64,109,35) .- 0.5) .* bestlogS

  end
end

# do the search in parallel
# syntax:  pmap((a1,a2)->f(a1,a2),{"a","b","c"},{2,1,3})
@time outarray = pmap((startθ,startη,startγ,startlogT,startlogS)->cenf.solveWithθ_f(startθ,startη,startγ,startlogT,startlogS,""),startθvec,startηvec,startγarr,startlogTarr,startlogSarr)

for i=1:CONST_SECONDTRIES
    ms_fvals[i] = outarray[i]["fval"]
    ms_θ[i] = outarray[i]["θ"]
    ms_β[i,:] = outarray[i]["β"]
    ms_η[i,:] = outarray[i]["η"]
    ms_γ[i,:] = outarray[i]["γ"]
    ms_logT[i,:] = outarray[i]["logT"]
    ms_logS[i,:] = outarray[i]["logS"]
    ms_r2[i] = outarray[i]["pseudor2"]
    try
      ms_varβ1[i] = outarray[i]["Var_β1"]
      ms_varβ2[i] = outarray[i]["Var_β2"]
    catch
      ms_varθ[i]= outarray[i]["Varθ"]
    end
    #CONST_DIM_ETA
    ms_varη1[i] = outarray[i]["Var_η1"]
    ms_varη2[i] = outarray[i]["Var_η2"]
    # ms_varη3[i] = outarray[i]["Var_η3"]
end
save("$(root)/output/ms_f_output_withvartheta_refined.jld2", "ms_fvals", ms_fvals, "ms_θ", ms_θ, "ms_β", ms_β, "ms_η", ms_η, "ms_γ", ms_γ, "ms_logT", ms_logT, "ms_logS", ms_logS)
